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Abstract 

A new approach to describe comminution processes in general ball mills as a 
macroscopic canonical ensemble is proposed. Using hamiltonian method, the 
model is able to take simultaneously into account the internal dynamics from 
mechanical motions inside the vial and external effects like electromagnetic 
and gravitational forces. Relevant physical observables are extracted using 
statistical mechanics approach through partition function at finite tempera- 
ture. The method enables numerical calculation using Monte Carlo technique 
to obtain, for instance particle number evolution in term of system tempera- 
ture. It is argued that the method is experimentally more verifiable than the 
conventional approaches based on geometrical displacements. An example of 
simulation for typical spex mill is also given. 

Key words: comminution, modeling, ball mill, hamiltonian, canonical 
ensemble 



1. Introduction 

The comminution processes in recent years attract the attention among 
scientists and engineers due to the increasing demand of ultrafine powders for 
nanotechnology applications in many areas. The demand then requires the 
improvement of comminution equipments like ball mills, roller mills and so 
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on. Unfortunately, the development of such comminution equipments always 
contains a lot of uncertainties due to a wide range of unknown parameters. 
These, in fact, lead to significant statistical errors. In order to overcome 
such problems, several models have been developed to quantitatively describe 
comminution process in various types of mills [1, 2, 3, 4, 5]. 

On the other hand, mathematical modeling and simulation may provide 
prior information and constraint to the unknown parameter ranges which 
should be useful to develop more optimized experimental strategy in com- 
minution processes. However, in most cases of mathematical models, the 
physical observables like grain-size etc are extracted from a set of equation 
of motions (EOMs). Such EOM's are considered to govern as complete as 
possible the dynamics of the system, from the mechanical motions to the 
evolution of grain-size distribution. This approach is obviously suffered from 
the nonlinearities of the equations under consideration, and then the require- 
ment of high computational power to solve them numerically. This fact often 
discourages a quantitative and deterministic approach for the simulation of 
such system. These nonlinear effects like chaotic behavior of the sphere 
motions within the mill encourages some works modeling the system using 
semi-empirical approaches [6]. However most of semi-empirical models re- 
quire a large number of experimental data based on prior observations [7] , or 
measured variables obtained from simulation results by other authors [8, 9]. 

More empirical approach is based on the physically realistic modeliza- 
tion of the ball mill system [10, 11]. This approach in general deals with 
three aspects : (1) evaluation of milling bodies dynamics and energetic in- 
puts transferred to powders; (2) description of the effects of such inputs on 
powders breakage; (3) description of powders evolution in terms of particle 
size distribution [12]. In a recent work [12], a comprehensive study on this 
line for the case of spex mixer / mill was performed by deploying the 3D 
simulation for milling bodies motion, and the population balance method to 
describe the granulometric evolution. Then, both are related through the 
energetic inputs in the population balance formalism which is obtained from 
the simulation of milling bodies motion. 

In this paper we propose a novel model and approach combining the 
deterministic approach for milling bodies motion, and the statistical approach 
to relate them with considerable macroscopic physical observables. This 
work is devoted to overcome the following problems in some conventional 
approaches : 
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• Experimentally it is almost impossible to trace the geometrical dis- 
placements of all matters in a vial with proper time resolution to verify 
the models which are based on the classical EOM. This problem is get- 
ting worse as one simulates a system with matters at few nanometers 
scale with comparable size of time-space resolution. 

• Solving a set of EOMs numerically, and then performing a simulation 
with high accuracy (enough time resolution) require huge efforts on 
both computing capacity and running time. 

• Taking into account the external circumstances around the vial like 
electromagnetic field and so forth. This might be interesting when 
one considers a comminution process which can reach the level of few 
nanometers. 

Therefore this work is intended to provide a tool for a nanometer system in 
a vial by developing direct relations between the vial internal dynamics with 
some external physical observables which should be more measurable. We 
should remark here that the vial internal dynamics is yet described empiri- 
cally using physical modelization approach. 

Further, rather solving the EOM's governing the whole dynamics, we use 
the hamiltonian approach to model all relevant interactions, and extract the 
physical observables through partition function by treating the system as a 
canonical ensemble at finite temperature. As a consequence, instead of ob- 
serving the geometrical evolution of matters in term of time in a ball mill, we 
can investigate the particle number distribution in term of temperature. This 
introduces a novel method relating the internal dynamics with the macro- 
scopic physical parameter like temperature, rather than time and geometrical 
displacements which are in most cases difficult to realize. The preliminary 
work but with incomplete hamiltonian on this matter has been reported in 
our previous work [13], however the present paper comprises more complete 
theoretical formalism including the external forces and the real simulation 
using Monte Carlo integration method. 

The paper is organized as follows. First, after this introduction we present 
the model and explain the underlying knowledge in detail. Before summa- 
rizing the results, numerical analysis and simulation for the case of ball mill 
with a structure similar to the well-known spex mixer / mill are discussed. 
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2. The model 



The whole system is modeled empirically using hamiltonian method. 
First we construct the total hamiltonian describing the dynamics in the ball 
mills. It is further followed by formulating the partition function and ex- 
tracting the relevant thermodynamics observables. 

2.1. The dynamics 

In our model, the dynamics of each 'matter' in the system, i.e. balls and 
powders inside the vial, are described by a hamiltonian H„^{f,t). The index 
m denotes the powder (p) or ball (b) and r = {x,y,z). The hamiltonian 
contains some terms representing all relevant interactions working on the 
matters inside the system as follow, 

Hm — Ho + Ki-m + Kn-v + Ki-m' + Kxt , (l) 

with V denotes the vial, while Hq is the free matter hamiltonian containing 
the kinetic term, 

1=1 

where is the matter number, and Pm are the matter mass and mo- 
mentum respectively. Throughout the paper we assume that the mass or size 
evolution of matters is uniform for the same matters. 

The matter self- interaction Vm_m, the matter-vial interaction Kn_v and 
the interactions between different matters may be induced by, for instance, 
impact (y™P) and Coulomb potentials, 

"m "m' /. ff ,).^ 

^) = - E E / d (e..O., n ■ (f^,) . . , (3) 

rim "m' 



1=1 j= 



1 \{fm)i - (fm'), 



with is the matter charge, while m, m' : v, p, b and n is the unit nor- 
mal vector. These potentials are considered describing the mechanical and 
static electrical properties of the matters. The first potentials should in fact 
represent the whole classical dynamics among the matters, i.e. the impact 
forces among balls and powders. This form will be clarified soon below. The 
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Coulomb potential disappears if the interacting matters have neutral charges. 
Also it works only in a short range of distance. Therefore, it should be negli- 
gible for some physics at nanometers scale as considered in the present case. 
The impact forces between the vial surface and balls or powders are treated 
in the same way using Eq. (3) by considering that the surface is formed by a 
set of fixed spheres [12]. This is inline with the simulation in the last section 
where the space displacement in the vial is at a distance comparable with 
the desired powder size, i.e. few tens nanometers at the largest. 

On the other hand, instead of Eq. (3) we can consider a simpler 'effective 
potential' hke the harmonic oscillator V^_^^, (r) = | /cmm' Ar^ to represent the 
whole mechanical properties in terms of 'effective coupling' /Cmm'- In this case, 
if rrira ^ 'iTT^m', the potential can be rewritten as Kn-m' = | i^m^ Ar^ since 
cUm = ^/ kj^m' / "tTT-m- Actually this is the case of ball-powder interactions since 
mb ^ rUp by the order of namely O(IO^). Nevertheless, in contrast with 
its simplicity, V°^^ absorbs the time dependencies and another interesting 
physical parameters characterizing the material properties like viscoelasticity. 
Young modulus etc. The time dependency is important to directly relate the 
system temperature with physical observables through finite temperature 
partition function as shown in the next subsection. Therefore, in this paper 
we take the impact potential to represent the mechanical properties in the 
system. 

The impact potential in Eq. (3) is induced by the impact force consists of 
normal and tangential components [12], F^,{f,t) — -^mm'(^'0 + -^mm'(^'0- 
The normal component may be written [14], 



2Tmm' /peff 3/2 3 /- 



n . (5) 



Here the first term is the elastic part based on the Hertz contact law, and 
the second term is the dissipative one that takes into account material vis- 
cosity. Tinni' is the Young modulus and Vnun' represents the Poisson ratio 
of the sphere material. The term R^^^ = {RmRm')/ {Rm + Rui') represents 
the effective radius, while ^mm' = Rm + Rm' — |^m — ^m'| is the displacement 
with i?m is the radius of interacting matter. A is a dissipative parameter 
[14, 15, 16], 



A 



1 3?7in' - Vi 



2 



m 



3 3?7m' + 



l-<m')(l-2^^: 



mm' 



J- mm' t^mm' 



(6) 
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The viscous constants r/m and rj^' relate the dissipative stress tensor to the 

deformation tensor [14, 15]. 

There are several proposed formulations for the tangential component 
F^^,{f,t). However it always assumes that the material slide upon each 

other in the case where the Coulomb condition ^ ^ram' — ^mm' holds, 
otherwise some viscous frictions occur [17]. In particular it follows that 
-^mm' ^) ^mm' ; "where the effective mass is m'^^, = m^m^i / {m^ + m^i ) 
[12]. Obviously, in our case with large mass discrepancy between powder and 
ball, i.e. mp/my, ~ which leads to ~ mp, the tangential impact force 
may be neglected for a good approximation. So, let us from now consider the 
normal component dominated impact force, that is F^,{f,t) ~ F^^,{f,t). 
This result simply yields the impact potential as written in Eq. (3) due to 
the Euler-Lagrange equation, 

F- ^ ^^^^ m 

since the dependency on matter velocity appears only in the tangential com- 
ponentF^^,(f,0 [12]. 

Beside the interactions among the matters itself, it is also possible to 
take into account the external potentials working on the whole system. For 
instance concerning the ball dynamics with relatively large ball size, the 
gravitational potential, 

rim 

VSr-rn^Gj2M, , (8) 

i=l 

might be important. Here, G is the gravitational constant. On the other 
hand, this should be less important for the powder dynamics due to its tiny 
size. 

On the other hand, we may also incorporate the effect of external electro- 
magnetic field surrounding the system to the charged matters. The potential 
is induced by the Lorentz force, F^^ — (E + v^n x -B), which leads to. 



EM 



rim 

= Qm5][0-(^'m)r^] , (9) 



and satisfies Eq. (7). (j) and A are the electromagnetic scalar and vector 
potentials related to the electric and magnetic fields by £^ = — V0 — dA/dt 
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and B = V X A. The inclusion of electromagnetic potential shifts the kinetic 
term in Eq. (2) as follow, 

Ho >Ho+EM=-^ ^\{Pm)i-QmA + Qm , (10) 

1=1 

Prom now, let us focus only on the dynamics of powders which is our 
main interest in the sense of comminution process. Prom Eqs. (1), (2), (3), 
(4) and (9), the total hamiltonian for the powder in our model is. 



1 I -2 

p 1=1 

-2 E E/ d(epp),,n.(F;7) 



-EEE/ d(U), n-(i^pZ),. , (11) 

m:b,v i=l j=l '^'^ 

for Qp 7^ 0. The last two potentials represent the total impact potential 
among powders; powders and vial; powders and balls respectively. Obviously 
we do not need to take into account the ball self-interaction V^^^b ^'^^ ball- 
vial interaction V^^^. This is actually the advantage of using hamiltonian 
method. 

2.2. Physical observables 

As mentioned briefly in introduction, the greatest advantage of deploying 
the hamiltonian method is one can extract some physical observables without 
solving the EOM's governing the system. Instead, in a canonical ensemble 
of matter m one can consider the partition function. 



/ ''-m pp 
Y\ dpi dfi exp - / dtHr, 
i=i L ^0 



(12) 



governed by a particular hamiltonian ifm- Here, /3 = l/lksT) with ks and 
T are the Boltzman constant and absolute temperature. Having partition 
function at hand, we can obtain some thermodynamics quantities in the 
system through relations namely, 

F^ = -^\nZ^, (13) 
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for free energy and, 



Prr. = 



' dV 



■ V ' 



(14) 



for pressure in a vial with volume V. 

In order to see the contributions from the interactions, it is more conve- 
nient to consider the normalized partition function. 



ZL = 



n 

i=i 



dpi dfi exp 
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dtH„ 



Zo 

KJrr 



n 

1=1 



dpi exp 



- dtHo^ 



(15) 



and further. 



PL^FL = 



InZn, 



(16) 



Performing the integral over time (t), we immediately obtain a tempera- 
ture dependent partition function, while the integrals over pi are the decou- 
pled gaussian integral which can be easily calculated. In the case of Eq. (11) 
it gives. 



2 mpTT 



/3 



np/2 



Y\ dfi exp 



i=l 



- / dtH' 



and. 



J L ^0 



(17) 



(18) 



where the interaction hamiltonian is. 



i(^i) = l 7 = 1 



EEE 

m:b,v i=l J=l 







(19) 



Obviously, only the scalar potential of external electromagnetic field con- 
tributes to the total energy of system under consideration. In other words, 

— * 

we can conclude here that in our model the magnetic field B does not infiu- 
ence the ball mill system, but the electric field E does. 
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Moreover, we can perform the integration over time (t) and to obtain 
further, 



4 - 



Y\ dfi exp 



T 



pp 



rim 



15 ^ 1 

m:b,v i=l j=l 



pm 



pm 



.5/2 
pm Vspin/ij 



(20) 



Prom this result, the thermodynamics observables are clearly not affected 
with the dissipative term, i.e. the second term in Eq. (5). After performing 
same integration we obtain. 



=1-/3^ In- (a^) , 



(21) 



respectively with. 



= 2 



Qp4>- 



rip nn 



15 n 



Up Up 

p i{^j)=i j=i 



T 



pp 



EEEr^^V^Kp".) 

•n:h.v i=l 1=1 P™ 



15 rir, 



PP 

5/2 
ij 



ReS ft ^5/2 

-"^pp vspp^j 



m:b,v i=l J=l 



(22) 



Eq. (21) provides a general behavior for temperature-dependent pressure in 
the model, while the geometrical structure and motion of vial is absorbed 
in the function T. Prom Eq. (21) clearly the physically meaningful regions 
are for < T < {2m-pT: ks)'^ and T > Tth- The later is equivalent to the 
condition, 

< ks Tth In (2 rup vr ks Tth) , (23) 

and Tth is always greater than (2 rup vr /cb)~^. The behavior of the temperature- 
dependent pressure in the model is depicted in Pig. 1. 



3. Simulation 

The simulation is done for a typical case of ball mill, that is spex mixer 
/ mill. However, the model under consideration can in principle be extended 
to deal with another types of ball mills by changing the coordinate system 
accordingly. 
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< r < (2 m KkJ 
T> 7"., 



T 



Figure 1: General contour of the normalized pressure as a function of system temperature. 



3.1. Coordinate system 

In the case of spex mill, the system contains a vial moving in a non- 
inertial system [10, 12]. Therefore, we can adopt the same coordinate system 
and schematization of vial surface as Figs. 1 and 2 in [12]. Since the system 
under consideration is moving on a non-inertial system f = {x,y,z), we 
should transform all coordinates in our formula to the coordinate system of 
inertial system {X, Y, Z). In the simulation, we make use of previous results 
on the roto-translation transformation [10, 12], 

X{t) — X cos9{t) cosa(t) +y cos 6{t) sin a{t) 

+z sin a{t) + L sine (t) , (24) 
Y{t) = -X sin a{t) + y cos a{t) , (25) 
Z{t) — —X sin^(t) cosQ;(t) —y sin9{t) sin a{t) 

+z cos a{t) + L cos 6(1) . (26) 

Here, L is the length of mechanical shaft-arm, 9 and a are the angles of 
rotation around the Y— and axis. Both angles can be written as follows, 

e = 00 sin{ujt + (f) , (27) 
a = ao sin{uj t + ip) , (28) 
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where 6q and ao are the angular momentum around the respected axis, co is 
the frequency and (p is the phase factor which depends on the initial condi- 
tions. 

In our simulation, for the sake of convenience the Cartesian coordinate is 
transformed into the cylindrical coordinate system, 

{x, y, z) — > (x, r cos 'd, r sin 'd) , (29) 

with 1? is the rotation angle around x in vial bases and r = -y/y^ + z^. 

As already mentioned above, we can apply the proposed model to any 
types of ball mill. This can be accomplished by replacing the coordinate 
system of ball mill under consideration like Eqs. (24)~(26) to the appropriate 
ones which represent its geometrical motion. 

3. 2. Technique 

The simulation within the present model and its underlying numerical 
calculation are done using Monte Carlo technique. This is the most appro- 
priate technique to deal with higher dimensional integral of many matters 
involved in the system like Eq. (22) [18]. 

We should note that the simulation here is performed to provide a com- 
plete picture on the model and its applications rather than showing a com- 
prehensive numerical simulation that is out of the scope of present paper. A 
comprehensive Monte Carlo integration in the present model requires proper 
resolution on 3-dimensional space which should be comparable with the pow- 
der size, and also enough time resolution within full running period to repre- 
sent the whole dynamics. Unfortunately, this kind of simulation is very time 
consuming, while the simulation itself is not the main interest in the present 
paper. As an illustration, for 100/xm powder size the appropriate resolution 
should be at least ~ 10^ x 10^ x 10^. The resolution must be increased ac- 
cordingly as the powder size is decreasing. Also, the present simulation is 
done for a particular time, that is equivalent to a particular static position 
of vial motion. Note that in the full simulation, the time evolution and the 
frequency of vial rotation are related each other. 

3.3. Results 

Now we are ready to perform a preliminary simulation for the powders 
with various strengths of interactions characterized by the defined parameters 
inside the potentials, and also in a circumstances with non-zero electric field 
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50 100 150 200 

Figure 2: as a function of the ratio of matter number, np/nb for silica, steel and nickel 
powders for R^/R^ = 50/3. 

— * 

E as well. The simulation is done for vial length / = 50 mm, vial radius 
Tv = 10 mm, shaft-arm length L = 200 mm and ball radius i?b = 5 mm. 

The simulation of J-" is performed for various values of the ratio of matter 
number {np/riy,) and matter size (i?b/-Rp)- Each case is also simulated for 
various materials of powders characterized with Young modulus (T) and 
Poisson ratio {v) of its sphere materials. The results are depicted in Figs. 2 
and 3. 

Figs. 2 and 3 show logarithmic scale of J-" for silica, steel and nickel 
powders. The error bars are coming from statistical errors due to the uncer- 
tainties of Young modulus and Poisson ratio of each material. The errors are 
significant for small ratio of matter number and large ratio of matter size. 
These facts are natural since as large as the matter number ratio, and also as 
small as matter size ratio would increase the probability of collisions between 
balls and powders. On the other hand, the increasing matter number and 
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10 20 30 40 50 10 20 30 40 



10 20 30 40 50 

R^/R 

Figure 3: J" as a function of the ratio of matter size, Rh/Rp for silica, steel and nickel 
powders for rip/nb = 30. 

size ratios indicate the on-going comminution processes inside the viaL 

Again, it should be remarked that the present simulation has lost the 
temperature evolution since it is done only for a single point of vial position. 
The temperature evolution will be recovered in a full simulation within the 
whole running period of ball mill. 

4. Summary 

We have proposed and discussed a novel model and approach for top-down 
mechanical comminution processes to produce nanomaterial using ball mill 
equipments. The study is focused on investigating the relevant potentials 
in the hamiltonian for a ball mill system and the formalism to extract rel- 
evant physical observables without tracking the geometrical displacements 
inside the vial. The governing hamiltonian is related to relevant physical 
observables through partition function of statistical mechanics approach. 
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Prom theoretical formalisms developed in this paper, we point out some 
interesting remarks which hold for any type (and geometrical motions) of 
ball mills, 

• Any ball mills should have the same temperature dependencies of its 
normalized pressure as shown in Eq. (21). Because the geometrical 
structure is absorbed in the auxiliary function T in Eq. (22). 

• The magnetic field does not affect the matter dynamics inside the vial, 
while the electric field does. 

• The contribution of dissipative term in the impact force is negligible 
for small ratio of mp/mb that is almost the case in all ball mills. This 
fact simplifies the whole analysis since the pre-knowledge of material 
viscosity is not needed anymore. 

• The model does not involve the powder number distribution, for in- 
stance using the population balance algorithm etc [19], in the formal- 
ism. The static powder number is represented by the number of in- 
teracting matters inside the vial, nm. Nevertheless, one can consider 
dynamic powder number distribution using any breakage functions like 
^raif) ~ ^Om(^o) exp(t/io) with n^^ is the initial powder number. 

According to these results, it can be argued that the model should eas- 
ily be confirmed by the experiments by verifying the above features. More 
importantly, the experimental setup is much simpler than the conventional 
approaches due to no requirement on tracking the geometrical displacements 
of involving matters inside the vial. 

We have also shown the simulation for particular case of spex mill us- 
ing Monte Carlo technique to calculate the complicated integrals in the for- 
malism. More detail and comprehensive simulations are still under progress. 
Such works involves some improvements. First of all, increasing the accuracy 
up to the realistic powder size at the order of ~ few tens nanometers. This 
means the Monte Carlo integration is done at the resolution of 10^ x 10^ x 10^. 
Secondly, performing complete temperature (time) evolution with appropri- 
ate accuracy within a realistic running period for various powder number 
distributions. According to the present preliminary result, the proper time 
resolution is ~ 0{t{/uj). 
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